Ensayo de Segmentación de imágenes RGB, basado en el ejercicio publicado por Ivan Lizarazo en formato de Notebook.
library ("rasterVis") |> suppressPackageStartupMessages()
library ('OpenImageR') |> suppressPackageStartupMessages()
library ('RImagePalette') |> suppressPackageStartupMessages()
library("raster") |> suppressPackageStartupMessages()
source(file = "indices_vegetacion_RGB.R")
##1. Leer la imagen para analizar
La imagen ha sido descargada del sitio web del Observatorio de la Tierra de la NASA . Ilustra el diverso paisaje agrÃcola en la parte occidental del estado de Minas Gerais en Brasil.
Se hace la asignación de la ruta para generalizar
# asignar la imagen con su ruta
# imagen <- "~/Pictures/balandra.jpg"
# imagen <- "~/giserias/conchalito-vuelo2/conchal_10Q.tif"
imagen <- "~/giserias/conchalito-vuelo2/RPubs_crops.jpg"
Se lee la imagen a analizar:
(crops <- stack(imagen)) # es un RasterStack con estructura para Ãndices
## class : RasterStack
## dimensions : 960, 1440, 1382400, 3 (nrow, ncol, ncell, nlayers)
## resolution : 1, 1 (x, y)
## extent : 0, 1440, 0, 960 (xmin, xmax, ymin, ymax)
## crs : NA
## names : RPubs_crops_1, RPubs_crops_2, RPubs_crops_3
Vemos la imagen cargada:
plotRGB(crops, 1, 2, 3) # visualiza un RasterStack
## Warning: [rast] unknown extent
## Warning: [rast] unknown extent
## Warning: [rast] unknown extent
Thanks to Matthias Forkel we can use a nice palette for vegetation indices:
colores <- function(n) {
.fun <- colorRampPalette(c("chocolate4", "orange", "yellow", "grey", "green", "green3", "darkgreen"))
col <- .fun(n)
return(col)
}
Let’s calculate the EG index:
# For the crops image red = 1, green=2, blue = 1.
(egvi <- EGVI(crops, 1, 2,3)) # los Ãndices regresan un RasterLayer
## Warning: [rast] unknown extent
## Warning: [rast] unknown extent
## Warning: [rast] unknown extent
## Warning: [rast] unknown extent
## Warning: [rast] unknown extent
## Warning: [rast] unknown extent
## class : RasterLayer
## dimensions : 960, 1440, 1382400 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 0, 1440, 0, 960 (xmin, xmax, ymin, ymax)
## crs : NA
## source : memory
## names : layer
## values : -0.3543046, 1.33945 (min, max)
Visualizamos:
cortes <- seq(-0.35, 1, by=0.1)
cols <- colores(length(cortes)-1)
# plot() despliega un RasterLayer
plot(egvi, col = cols, breaks=cortes, main = "EG Vegetation Index")
Calculamos el Ãndice NGRDI:
# For the crops image red = 1, green=2, blue = 1.
(ngvi <- NGRDI(crops, 1, 2,3))
## Warning: [rast] unknown extent
## Warning: [rast] unknown extent
## Warning: [rast] unknown extent
## Warning: [rast] unknown extent
## class : RasterLayer
## dimensions : 960, 1440, 1382400 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 0, 1440, 0, 960 (xmin, xmax, ymin, ymax)
## crs : NA
## source : memory
## names : layer
## values : -0.3809524, 1 (min, max)
Visualizamos:
cortes <- seq(-0.35, 1.0, by=0.1)
cols <- colores(length(cortes)-1)
##plot(ndvitrend, 2, col=cols, breaks=classbreaks)
plot(ngvi, col = cols, breaks=cortes, main = "NGRDI Vegetation Index")
Leamos la imagen usando el paquete OpenImageR:
#RGB.Color <- RGB.Color[, , 1:3] # si lee cuatro capas en lugar de tres
RGB.Color <- readImage(imagen) # Large array
# Nota: las funciones readJPGE, readPNG leen a matriz la imagen, tal y como lo hace
# {OpenImageR}::readImage() lee a matriz o array los formatos .png, .jpeg, .jpg, .tif .tiff
Usemos ahora SLIC para la segmentación de imágenes. Para obtener una explicación detallada del código, consultar la viñeta del paquete.
# superpixel recibe un Large array o imagen en forma de matriz con tres capas de fondo
Region.slic = superpixels(input_image = RGB.Color, method = "slic", superpixel = 80,
compactness = 30, return_slic_data = TRUE,
return_labels = TRUE, write_slic = "",
verbose = FALSE)
## Warning in interface_superpixels(input_image, method, superpixel, compactness,
## : The input data has values between 0.000000 and 1.000000. The image-data will
## be multiplied by the value: 255!
imageShow(Region.slic$slic_data)
str(RGB.Color)
## num [1:960, 1:1440, 1:3] 0.886 0.855 0.843 0.82 0.749 ...
Let’s get a matrix with NGRDI values:
# ngvi
# egvi
# egvi@data@values
# ngvi@data@values
ngvi.mat <- matrix(ngvi@data@values,
nrow = ngvi@nrows,
ncol = ngvi@ncols, byrow = TRUE)
# min(ngvi.mat)
# max(ngvi.mat)
La función imageShow() trabaja con datos que están en el rango de ocho bits 0 – 255 o en el rango [0,1] (es decir, el rango de x entre 0 y 1 inclusive). Sin embargo, no funciona con valores NGRVI si estos valores son negativos. Por lo tanto, escalaremos los valores de NGRDI a [0,1].
#rescale between 0 and 1; ver la funcion {OpenImageR}::NormalizeObject() que se usa mas abajo
normalize <- function(x) {
min <- min(x)
max <- max(x)
return(1* (x - min) / (max - min))
}
ngvi.mat_norm <- normalize(ngvi.mat)
imageShow(ngvi.mat_norm)
ngvi.data <- RGB.Color
# ngvi.data[,,1] <- ngvi.mat_norm
# ngvi.data[,,2] <- ngvi.mat_norm
# ngvi.data[,,3] <- ngvi.mat_norm
En la siguiente sección describiremos cómo crear una segmentación de imágenes de los datos del NGRDI y cómo utilizar el análisis de conglomerados para crear una clasificación del terreno
Aquà se muestra una aplicación de la función superpixels() para la clasificación de la cobertura terrestre de los datos NGRDI.
ngvi.80 = superpixels(input_image = ngvi.data,
method = "slic", superpixel = 80,
compactness = 30, return_slic_data = TRUE,
return_labels = TRUE, write_slic = "",
verbose = TRUE)
## Warning in interface_superpixels(input_image, method, superpixel, compactness,
## : The input data has values between 0.000000 and 1.000000. The image-data will
## be multiplied by the value: 255!
## The input image has the following dimensions: 960 1440 3
## The 'slic' method will be utilized!
## The output image has the following dimensions: 960 1440 3
imageShow(ngvi.80$slic_data)
Aquà la estructura del objeto ngvi.80.
str(ngvi.80)
## List of 2
## $ slic_data: num [1:960, 1:1440, 1:3] 226 218 215 209 191 169 151 136 130 127 ...
## $ labels : num [1:960, 1:1440] 0 0 0 0 0 0 0 0 0 0 ...
Es una lista con dos elementos, slic_data, una matriz tridimensional de datos de color de pÃxeles (no normalizados) y etiquetas, una matriz cuyos elementos corresponden a los pÃxeles de la imagen. El nombre del segundo elemento sugiere que puede contener valores que identifican el superpÃxel al que pertenece cada pÃxel. Veamos si esto es cierto.
sort(unique(as.vector(ngvi.80$labels)))
## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## [26] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## [51] 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73
Hay 72 etiquetas únicas. Aunque la llamada a superpixels() especificó 80 superpÃxeles, la función generó 72. Podemos ver qué pÃxeles tienen el valor de etiqueta 0 estableciendo los valores de todos los demás pÃxeles en [255, 255, 255], lo que se trazará como blanco.
R0 <- ngvi.80
for (i in 1:nrow(R0$label))
for (j in 1:ncol(R0$label))
if (R0$label[i,j] != 0)
R0$slic_data[i,j,] <- c(255,255,255)
imageShow(R0$slic_data)
La operación aÃsla el superpÃxel en la esquina superior derecha de la imagen, junto con la porción correspondiente del lÃmite. Podemos utilizar fácilmente este enfoque para determinar qué valor de ngvi.80$label corresponde a qué superpÃxel. Ahora tratemos con el lÃmite. Una pequeña exploración del objeto ngvi.80 sugiere que los pÃxeles en el lÃmite tienen los tres componentes iguales a cero. Aislamos y trazamos todos esos pÃxeles coloreando todos los demás pÃxeles de blanco.
Bdry <- ngvi.80
for (i in 1:nrow(Bdry$label))
for (j in 1:ncol(Bdry$label))
if (!(Bdry$slic_data[i,j,1] == 0 &
Bdry$slic_data[i,j,2] == 0 &
Bdry$slic_data[i,j,3] == 0))
Bdry$slic_data[i,j,] <- c(255,255,255)
Bdry.norm <- NormalizeObject(Bdry$slic_data)
imageShow(Bdry$slic_data)
La figura muestra que efectivamente se han identificado los pÃxeles lÃmite. Hay que tener en cuenta que la función imageShow() muestra estos pÃxeles en blanco con un borde negro, en lugar de negro puro.
Habiendo realizado un análisis preliminar, podemos organizar nuestro proceso de segmentación en dos pasos.
El primer paso será reemplazar cada uno de los superpÃxeles generados por la función superpixels() de OpenImageR por uno en el que cada pÃxel tenga el mismo valor, correspondiente a una medida de tendencia central (p. ej., la media, mediana o moda) del superpÃxel original.
El segundo paso será utilizar el procedimiento de agrupamiento no supervisado de K-medias para organizar los superpÃxeles del paso 1 en un conjunto de grupos y darle a cada grupo un valor correspondiente a una medida de tendencia central del grupo.
Utilizaremos la función make.segments() para realizar la segmentación. El primer argumento de make.segments() es el objeto de superpÃxeles y el segundo es la medida funcional de tendencia central. Aunque en este caso cada uno de los tres colores del objeto ngvi.80 tiene los mismos valores, esto puede no ser cierto para todas las aplicaciones, por lo que la función analiza cada color por separado.
La función, escrita por Richard Plant, es la siguiente:
# Identify a measure of central tendency of each superpixel
make.segments <- function(x, ftn){
# The argument ftn is any functional measure of central tendency
z <- x
# For each identified superpixel, compute measure of central tendency
for (k in unique(as.vector(x$labels))){
# Identify members of the superpixel having the given label
in.super <- matrix(0, nrow(x$label), ncol(x$label))
for (i in 1:nrow(x$label))
for (j in 1:ncol(x$label))
if (x$label[i,j] == k)
in.super[i,j] <- 1
#Identify the boundary cells as having all values 0
on.bound <- matrix(0, nrow(x$label), ncol(x$label))
for (i in 1:nrow(x$label))
for (j in 1:ncol(x$label))
if (in.super[i,j] == 1){
if (x$slic_data[i,j,1] == 0 & x$slic_data[i,j,2] == 0
& x$slic_data[i,j,3] == 0)
on.bound[i,j] <- 1
}
#Identify the superpixel cells not on the boundary
sup.data <- matrix(0, nrow(x$label), ncol(x$label))
for (i in 1:nrow(x$label))
for (j in 1:ncol(x$label))
if (in.super[i,j] == 1 & on.bound[i,j] == 0)
sup.data[i,j] <- 1
# Compute the measure of central tendency of the cells in R, G, B
for (n in 1:3){
# Create a matrix M of the same size as the matrix of superpixel values
M <- matrix(0, dim(x$slic_data)[1], dim(x$slic_data)[2])
for (i in 1:nrow(x$label))
for (j in 1:ncol(x$label))
# Assign to M the values in the superpixel
if (sup.data[i,j] == 1) M[i,j] <- x$slic_data[i,j,n]
if (length(M[which(M > 0 & M < 255)]) > 0)
# Compute the measure of central tendency
ftn.n <- round(ftn(M[which(M > 0 & M < 255)]), 0)
else
ftn.n <- 0
for (i in 1:nrow(x$label))
for (j in 1:ncol(x$label))
if (in.super[i,j] == 1) z$slic_data[i,j,n] <- ftn.n
}
}
return(z)
}
Aquà está la aplicación de esa función al objeto ngvi.80 con el segundo argumento configurado como mean.
Podemos ver el resultado.
ngvi.means <- make.segments(ngvi.80, mean)
imageShow(ngvi.means$slic_data)
El siguiente paso es desarrollar grupos que representen tipos de cobertura terrestre identificables. En un proyecto real, el procedimiento serÃa recopilar un conjunto de datos reales del sitio, pero esa opción no está disponible para nosotros. En su lugar, trabajaremos con la reproducción en color real de la foto del ejemplo.
La cobertura del suelo se puede subdividir utilizando K-medias en cualquier número de tipos. Tenga en cuenta que no se ha interpretado aquà la relación entre los grupos de K-medias y las clases de cobertura del suelo existentes en el sitio del ejemplo. Sin embargo, en una aplicación real, esta tarea es imprescindible.
set.seed(123)
# OJO: en la función kmeans se establece el número de clases subsecuentes
nclass <- 10
ngvi.clus <- kmeans(as.vector(ngvi.means$slic_data[,,1]), nclass)
vege.class <- matrix(ngvi.clus$cluster,
nrow = ngvi@nrows,
ncol = ngvi@ncols, byrow = FALSE)
class.ras <- raster(vege.class,
xmn = crops@extent@xmin,
xmx = crops@extent@xmax,
ymn = crops@extent@ymin,
ymx = crops@extent@ymax,
crs = crs(crops))
A continuación podemos usar la función {raster} ratify() para asignar niveles de factores descriptivos a los grupos o clases definidas.
class.ras <- ratify(class.ras)
rat.class <- levels(class.ras)[[1]]
# eventualmente, se asignan etiquetas para cada una de las clases del resultado;
# en el caso del ejemplo, se aplica:
# rat.class$landcover <- c("LC 10", "LC 09", "LC 08", "LC 07", "LC 06", "LC 05", "LC 04", "LC 03", "LC 02", "LC 01")
rat.class$landcover <- paste("Clase", nclass:1)
#Crear una palette con el número de colores correspondientes al resultado en kmeans
# se lee de nuevo la imagen original como raster comun en Large array x, y, 1:3
# ncrops <- jpeg::readJPEG("~/Pictures/balandra.jpg")
ncrops <- readImage(imagen) # para no invocar otra librerÃa
crops.pal <- image_palette(ncrops, n=nclass)
# Nota: las funciones readJPGE, readPNG leen a matriz la imagen, tal y como lo hace
# {OpenImageR}::readImage() lee a matriz o array los formatos .png, .jpeg, .jpg, .tif .tiff
levels(class.ras) <- rat.class
levelplot(class.ras, margin=FALSE,
col.regions= crops.pal[1:nclass],
# se pueden reasignar los colores asà como agrupar asignando el mismo a diferentes clases
# col.regions= crops.pal[c(5,7,9,4,8,10,2,4,1,6)],
main = "Land Cover Types in crops")